if(!require(OpenStreetMap)){install.packages("OpenStreetMap"); 
  library(OpenStreetMap)}

## Loading required package: OpenStreetMap

if(!require(ggplot2)){install.packages("ggplot2"); library(ggplot2)}

## Loading required package: ggplot2

* 웹에서 제공하는 지도정보를 사용할 수 있는 패키지
* `OpenStreetMap`을 이용하기 위한 패키지 호출

map = OpenStreetMap::openmap(upperLeft = c(43, 119), 
                             lowerRight = c(33, 134), type = 'bing')
plot(map)

    * upperLeft에는 지도 왼쪽 위 모서리의 위경도 좌표를 입력
    * lowerRight에는 지도 오른쪽 아래 모서리의 위경도 좌표를 입력
    * ggplot2 패키지의 autoplot 함수로 openmap 을 그릴 수 있음

nm = c("osm", "mapbox", "stamen-toner", 
       "stamen-watercolor", "esri", "esri-topo", 
       "nps", "apple-iphoto", "osm-public-transport")
par(mfrow=c(3,3),  mar=c(0, 0, 0, 0), oma=c(0, 0, 0, 0))

for(i in 1:length(nm)){
  map <- openmap(c(43,119),
                 c(33,134),
                 minNumTiles = 3,
                 type = nm[i])
  plot(map, xlab = paste(nm[i]))
}
par(mfrow = c(1, 1))

* openmap 함수의 type의 종류는 bing 외에도 많이 있다.
    + nm 벡터에 웹지도의 소스에 대한 이름이 들어 있다.

• mar, oma option : 그래프의 내부, 외부 마진 공간 조절해주는 옵션

map1 <- openmap(c(43.46,119.94),
               c(33.22,133.98))
plot(map1) 
abline(h = 38, col = 'blue')

str(map1)

## List of 2
##  $ tiles:List of 1
##   ..$ :List of 5
##   .. ..$ colorData : chr [1:378378] "#FAF9F6" "#585857" "#FAF9F6" "#FAF9F6" ...
##   .. ..$ bbox      :List of 2
##   .. .. ..$ p1: num [1:2] 13351660 5382253
##   .. .. ..$ p2: num [1:2] 14914585 3924542
##   .. ..$ projection:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. ..@ projargs: chr "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
##   .. ..$ xres      : int 594
##   .. ..$ yres      : int 637
##   .. ..- attr(*, "class")= chr "osmtile"
##  $ bbox :List of 2
##   ..$ p1: num [1:2] 13351660 5382253
##   ..$ p2: num [1:2] 14914585 3924542
##  - attr(*, "zoom")= int 6
##  - attr(*, "class")= chr "OpenStreetMap"

map1$tiles[[1]]$projection

## CRS arguments:
##  +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0
## +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

• Mercator(EPSG:3857)로 위경도 좌표계가 아님

if(!require(sp)){install.packages("sp"); library(sp)}

map_p <- openproj(map1, projection = CRS("+proj=longlat"))
str(map_p)

## List of 2
##  $ tiles:List of 1
##   ..$ :List of 5
##   .. ..$ colorData : chr [1:393336] NA NA NA NA ...
##   .. ..$ bbox      :List of 2
##   .. .. ..$ p1: num [1:2] 119.8 43.6
##   .. .. ..$ p2: num [1:2] 134.1 33.1
##   .. ..$ projection:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. .. ..@ projargs: chr "+proj=longlat +ellps=WGS84"
##   .. ..$ xres      : int 607
##   .. ..$ yres      : int 648
##   .. ..- attr(*, "class")= chr "osmtile"
##  $ bbox :List of 2
##   ..$ p1: num [1:2] 119.8 43.6
##   ..$ p2: num [1:2] 134.1 33.1
##  - attr(*, "zoom")= int 6
##  - attr(*, "class")= chr "OpenStreetMap"

plot(map_p)
abline(h = 38, col = 'blue')

• 우리나라 지역에 잘 맞는 위경도 변환

map_p <- openproj(map1, projection = 
                    CRS("+proj=utm +zone=52N + datum=WGS84"))
plot(map_p)
abline(h = 38, col = 'blue')

• 좌표계 변환 과정이 필요

좌표계변환 (sp 라이브러리의 활용)

• 좌표 데이터프레임 생성
• sp 클래스로 변환
• projection 설정
• 좌표계 변환
• 지도위에 매핑

• 좌표 데이터프레임 생성 및 sp 클래스 변환

if(!require(sp)){install.packages("sp"); library(sp)}

a  <-data.frame(lon =  seq(100,140,by = 0.1),
           lat =  38)
sp::coordinates(a) = ~ lon + lat
str(a)

## Formal class 'SpatialPoints' [package "sp"] with 3 slots
##   ..@ coords     : num [1:401, 1:2] 100 100 100 100 100 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:401] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] 100 38 140 38
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

• projection 설정

sp::proj4string(a) = "+proj=longlat"
#a@proj4string  = CRS("+proj=longlat")
str(a)

## Formal class 'SpatialPoints' [package "sp"] with 3 slots
##   ..@ coords     : num [1:401, 1:2] 100 100 100 100 100 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:401] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] 100 38 140 38
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=longlat +ellps=WGS84"

• 좌표계 변환 (sp::CRS)

a_tf = spTransform(a,  CRS("+proj=utm +zone=52N + datum=WGS84"))
str(a_tf)

## Formal class 'SpatialPoints' [package "sp"] with 3 slots
##   ..@ coords     : num [1:401, 1:2] -2069280 -2060288 -2051296 -2042306 -2033316 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:401] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] -2069280 4205815 1467197 4626502
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=utm +zone=52N +ellps=WGS84"

• 매핑

plot(map_p)
points(a_tf@coords[,1], a_tf@coords[,2], type = 'l', col = 'blue')

* `openmap`으로 그린 지도위에 새로운 시각정보를 그리기 위해서 
`openmap`이 생성한 좌표계를 알아야 한다. 
    + 위경도 
    + 구면좌표의 지도위로의 정사형 방법

* 새로운 시각정보를 표시하는 위치는 이미 그려진 지도의 좌표계를 따른다.

* 좌표계를 변환하는 라이브러리인 `sp`라이브러리를 사용하자.

if(!require(sp)){install.packages("sp"); library(sp)}
if(!require(mapplots)){install.packages("mapplots"); library(mapplots)}

map = openmap(upperLeft = c(43, 119),lowerRight = c(33, 134))
seoul_loc = geocode('Seoul')
coordinates(seoul_loc) = ~lon + lat
proj4string(seoul_loc) = "+proj=longlat +datum=WGS84" 
seoul_loc_Tf = spTransform(seoul_loc,
                  CRS(as.character(map$tiles[[1]]$projection)))
plot(map)
add.pie(z=1:2,labels = c('a','b'),
        x = seoul_loc_Tf@coords[1],
        y = seoul_loc_Tf@coords[2], radius = 100000)

* `geocode` 함수를 이용하여 seoul의 위경도 좌표 호출
* `coordinates(seoul_loc) = ~lon + lat`를 이용하여 `sp`라이브러리에서 처리가능한
데이터 형으로 변환
* `sp` 패키지의 `spTransform()` 을 이용, openmap의 좌표계로 위경도 좌표를 변환
    + `map$tiles[[1]]$projection`에 이미 읽어들인 openStreetmap의 좌표정보가 있음.
    +  `CRS`  함수를 이용하여  `CRSobj` 데이터 형으로 변환
    

if(!require(ggmap)){install.packages("ggmap"); library(ggmap)}

* ggmap을 이용하기 위한 패키지 호출
* 휴스턴의 범죄 데이터의 시각화 작업

data(crime)
head(crime, 2)

##                      time     date hour premise offense  beat     block
## 82729 2010-01-01 15:00:00 1/1/2010    0     18A  murder 15E30 9600-9699
## 82730 2010-01-01 15:00:00 1/1/2010    0     13R robbery 13D10 4700-4799
##          street type suffix number   month    day                 location
## 82729   marlive   ln      -      1 january friday    apartment parking lot
## 82730 telephone   rd      -      1 january friday road / street / sidewalk
##                 address       lon      lat
## 82729   9650 marlive ln -95.43739 29.67790
## 82730 4750 telephone rd -95.29888 29.69171

* crime 데이터에는 범죄 발생 시기, 위치, 종류 등이 저장

violent_crimes = subset(crime,
                   offense != "auto theft" & 
                   offense != "theft" & 
                   offense != "burglary")
violent_crimes$offense <- factor(violent_crimes$offense,
      levels = c("robbery", "aggravated assault", "rape", "murder"))
violent_crimes = subset(violent_crimes,
                      -95.39681 <= lon & lon <= -95.34188 &
                      29.73631 <= lat & lat <=  29.78400)

* 범죄 데이터를 필요한 것만 선택해서 네가지의 범죄 종류를 factor 형태로 저장
    + `levels`:범죄의 중대성 순서대로 level을 정했음
* Houston의 downtown에서 발생한 범죄에 한정

HoustonMap = qmap("houston", zoom = 14,
                color = "bw", legend = "topleft")
HoustonMap + geom_point(aes(x = lon, y = lat,
                            colour = offense, size = offense),
                            data = violent_crimes)

* `ggmap` 패키지의 `qmap` 함수를 이용하여 Houston의 지도 저장
    + `zoom`: 확대정도
    + `color`: 지도의 종류 

* `geom_point`를 이용해서 저장된 지도 위에 `violent_crimes`를 점으로 시각화
    + `colour = offense`: 색깔을 이용해서 범죄의 종류를 구분
    + `colour = offense`: 범죄의 종류를 점의 크기에 적용
       중대한 범죄의 레벨값을 높게 정했음

* Houston 지도 위의 범죄 데이터를 점으로 표기 <br>


* geom_density2d 함수를 이용하여 범죄 밀도 시각화
    + `size`: 선의 굵기
    + `bins`: 축방향의 bin 개수

HoustonMap +
    geom_point(aes(x = lon, y = lat,
                   colour = offense, size = offense),
               data = violent_crimes) +
    geom_density2d(aes(x = lon, y = lat), size = 0.2 , bins = 4, 
                   data = violent_crimes)

## Crimes map

* `stat_density2d` 을 이용한 density plot의 색깔 칠하기 

    load('airport.Rdata')
    head(airport_krjp)

##              city country iata      lat      lon Freq
## 1           Tokyo   Japan  NRT 35.76472 140.3864  100
## 2       Matsumoto   Japan  MMJ 36.16676 137.9227    7
## 3         Ibaraki   Japan  IBR 36.18108 140.4154    4
## 4         Iwojima   Japan  IWO 24.78400 141.3227   NA
## 5 Nanki-shirahama   Japan  SHM 33.66222 135.3644    2
## 6         Obihiro   Japan  OBO 42.73333 143.2172    6

* airport.Rdata는 `airport_krjp`, `link_krjp` 데이터를 포함
* airport_krjp는 공항의 데이터

head(link_krjp)

##   iata      lat      lon group
## 1  CJJ 36.71660 127.4991     1
## 2  CJU 33.51131 126.4930     2
## 3  CJU 33.51131 126.4930     3
## 4  CJU 33.51131 126.4930     4
## 5  FUK 33.58594 130.4507     5
## 6  GMP 37.55831 126.7906     6

* link_krjp는 비행의 출발지와 도착지의 데이터
* link_krjp 데이터에서 하나의 비행경로는 하나의 group으로 표기

map = ggmap(get_googlemap(center = c(lon=134, lat=36),
                          zoom = 5, maptype='roadmap', color='bw'))
map + geom_line(data=link_krjp,aes(x=lon,y=lat,group=group), 
                col='grey10',alpha=0.05) + 
  geom_point(data=airport_krjp[complete.cases(airport_krjp),],
             aes(x=lon,y=lat, size=Freq), colour='black',alpha=0.5) + 
  scale_size(range=c(0,15))

* google map에서 지도를 가져와서 map 변수에 저장
    + `google map`은 기본적으로 위경도 좌표계를 사용하므로 좌표변환이 필요없음
* 지도에 `geom_line` 함수를 이용하여 비행 경로를 시각화
* `geom_point` 함수를 이용하여 점의 크기에 따라 공항 이용 빈도를 시각화

* 지도에 몇 개의 점에서 관측치를 표시

* 관측되지 않은 지점의 예측값을 표시하고 싶은 경우
    + 공간적 내삽 (spatial interpolation)
    
* Kriging 을 이용한 서울시 미세먼지 농도의 시각화 작업

if (!require(sp)) {install.packages('sp'); library(sp)}
if (!require(gstat)) {install.packages('gstat'); library(gstat)}
if (!require(automap)) {install.packages('automap'); library(automap)}
if (!require(rgdal)) {install.packages('rgdal'); library(rgdal)}
if (!require(e1071)) {install.packages('e1071'); library(e1071)}
if (!require(dplyr)) {install.packages('dplyr'); library(dplyr)}
if (!require(lattice)) {install.packages('lattice'); library(lattice)}
if (!require(viridis)) {install.packages('viridis'); library(viridis)}

데이터 정제 및 시각화를 위한 library 호출

seoul032823 <- read.csv ("seoul032823.csv")
head(seoul032823)

##     X.1     X     ID       time PM10      LAT      LON      day
## 1   671   671 111121 2012032823  127 37.56464 126.9760 20120328
## 2  2879  2879 111123 2012032823  121 37.57203 127.0050 20120328
## 3  5087  5087 111131 2012032823  127 37.54031 127.0051 20120328
## 4  7295  7295 111141 2012032823  143 37.54464 127.0957 20120328
## 5  9503  9503 111142 2012032823  140 37.54311 127.0411 20120328
## 6 11711 11711 111151 2012032823  128 37.58495 127.0943 20120328

* 서울시에서 제공하는 미세먼지 데이터 읽기

skorea <- raster::getData(name ="GADM", country= "KOR", level=2)
# skorea <- readRDS("KOR_adm2.rds")
head(skorea,2)

## class       : SpatialPolygonsDataFrame 
## features    : 2 
## extent      : 128.9957, 129.0844, 35.12158, 35.28471  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 15
## names       : OBJECTID, ID_0, ISO,      NAME_0, ID_1, NAME_1, ID_2,   NAME_2, HASC_2, CCN_2, CCA_2, TYPE_2, ENGTYPE_2,          NL_NAME_2, VARNAME_2 
## min values  :        1,  213, KOR, South Korea,    1,  Busan,    1,      Buk,       ,    NA,      ,     Gu,  District, 부산진구| 釜山鎭區,           
## max values  :        2,  213, KOR, South Korea,    1,  Busan,    2, Busanjin,       ,    NA,      ,     Gu,  District,         북구| 北區,

• GADM에서 데이터 가져오기

데이터 출처 : Global Administrative Areas

class(skorea)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

head(skorea@polygons[[1]]@Polygons[[1]]@coords, 3)

##          [,1]     [,2]
## [1,] 128.9957 35.17683
## [2,] 128.9983 35.19250
## [3,] 129.0003 35.20673

* GADM 데이터 구조 파악
    +데이터의 구조가 'sp' 패키지에서 제공하는 S4-class 구조임을 확인
    +head 함수로 확인할 수 없었던 좌표계를 @ 을 이용하여 확인

if (!require(broom)) {install.packages('broom'); library(broom)}

## Loading required package: broom

skorea <- broom::tidy(skorea)

## Regions defined for each Polygons

class(skorea)

## [1] "data.frame"

head(skorea,3)

##       long      lat order  hole piece group id
## 1 128.9957 35.17683     1 FALSE     1   1.1  1
## 2 128.9983 35.19250     2 FALSE     1   1.1  1
## 3 129.0003 35.20673     3 FALSE     1   1.1  1

* 데이터의 구조를 data.frame 형태로 변경

ggplot() + geom_map(data= skorea, map= skorea,
  aes(map_id=id,group=group),fill=NA, colour="black") +
  geom_point(data=seoul032823, aes(LON, LAT, col = PM10),alpha=0.7) +
  labs(title= "PM10 Concentration in Seoul Area at South Korea",
       x="Longitude", y= "Latitude", size="PM10(microgm/m3)")

* geom_map 함수를 이용하여, GADM 데이터를 지도위에 plotting<br>
* geom_point 함수를 이용하여, 미세먼지 데이터를 점으로 표기<br>
* 보다 나은 시각화를 위해 점과 점 사이에 미세먼지 값의 내삽

class(seoul032823)

## [1] "data.frame"

coordinates(seoul032823) <- ~LON+LAT
class(seoul032823)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

 * 데이터 구조 변경
    +Kriging 방법 적요을 위해 데이터를 'sp dataframe' 형태로 변경
    + 위경도 정보는 ``seoul032823@coords``에 저장되어 있음

LON.range <- c(126.5, 127.5)
LAT.range <- c(37, 38)
seoul032823.grid <- expand.grid(
  LON = seq(from = LON.range[1], to = LON.range[2], by = 0.01),
  LAT = seq(from = LAT.range[1], to = LAT.range[2], by = 0.01))

* Kriging을 위한 grid 생성
    + 새로운 좌표계의 범위를 계산 
    + 좌표계의 범위를 이용하여 grid 생성

plot(seoul032823.grid)
points(seoul032823, pch= 16,col="red")

* 참고: grid 위에서 데이터 시각화

coordinates(seoul032823.grid)<- ~LON+LAT ## sp class
gridded(seoul032823.grid)<- T
plot(seoul032823.grid)
points(seoul032823, pch= 16,col="red")

    * grid 데이터의 구조 변경 및 격자로 변경후 시각화

seoul032823_OK <- autoKrige(formula = PM10~1,
                            input_data = seoul032823,
                            new_data = seoul032823.grid )

## [using ordinary kriging]

* `autoKrige` 함수를 통해 만들어직 격자데이터 `seoul032823.grid` 위에서
` seoul032823` 의 데이터를 내삽함

* 결과는 `seoul032823_OK`에 저장됨

head(seoul032823_OK$krige_output@coords, 2)

##      LON LAT
## 1 126.50  37
## 2 126.51  37

head(seoul032823_OK$krige_output@data$var1.pred,2)

## [1] 169.6544 169.6502

* 내삽 정보는위의 객체에 저장

myPoints <- data.frame(seoul032823)
myKorea <- data.frame(skorea)
myKrige <- data.frame(seoul032823_OK$krige_output@coords, 
              pred = seoul032823_OK$krige_output@data$var1.pred)   

* ggplot을 사용하기 위해 데이터 형식을 data.frame으로 변형

if(!require(viridis)){install.packages("viridis"); library(viridis)}

ggplot()+ theme_minimal() +
  geom_tile(data = myKrige, aes(x= LON, y= LAT, fill = pred)) +
  geom_map(data= myKorea, map= myKorea, aes(map_id=id,group=group),
           fill=NA, colour="black") +
  coord_cartesian(xlim= LON.range ,ylim= LAT.range) +
  labs(title= "PM10 Concentration in Seoul Area at South Korea",
       x="Longitude", y= "Latitude")+
  theme(title= element_text(hjust = 0.5,vjust = 1,face= c("bold")))+
  scale_fill_viridis(option="magma")